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Abstract 

In this letter we discuss the vahdity of the ergodicity hypothesis in theories of violent relaxation 
in long-range interacting systems. We base our reasoning on the Hamiltonian Mean Field model 
and show that the life-time of quasi-stationary states resulting from the violent relaxation does 
not allow the system to reach a complete mixed state. We also discuss the applicability of a 
generalization of the central limit theorem. In this context, we show that no attractor exists in 
distribution space for the sum of velocities of a particle other than the Gaussian distribution. The 
long-range nature of the interaction leads in fact to a new instance of sluggish convergence to a 
Gaussian distribution. 
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I. INTRODUCTION 



Physical systems with long-range interaction has been widely studied with results re- 
ported in the literature, and properties that are markedly different from short range in- 
teracting systems, e.g, non-Gaussian quasi- stationary states, negative heat-capacity, tem- 
perature jumps at critical points, anomalous diffusion and ensemble inequivalence 
Examples include self-gravitating systems, vortices in two-dimensional fluids, dipolar inter- 
actions, fractures in solids and simplified models as the Hamiltonian Mean Field (HMF), 
free electron laser and the plasma single wave models . The nature of the 

non-Gaussian quasi-stationary states arising from the dynamics can be interpreted as stable 
stationary states of the mean-field Vlasov equation which describes the system in the limit 

TIQ 

— oo, where N is the number of particles [8|, |9|. This equivalence was explicitly verified 
by comparing results from direct simulations of a hamiltonian system and the numerical 
solutions of the associated Vlasov equation For longer times correlations build between 
particles and the Vlasov equation must be replaced by a proper kinetic equation UJ, Il2 |. 
The initial stage of time evolution is called violent relaxation and corresponds to a very 
rapid evolution from the initial condition into a quasi-stationary state. Since the seminal 



paper by Lynden-Bell [13| , many attempts were made to formulate the statistical mechanics 



of violent relaxation (see for instance References [14[] and 15| and references therein), but 
a satisfactory theory for predicting the outcome of the violent relaxation process is still 
lacking. 

In the last decade many papers tried to identify the outcome states of a violent relaxation 
with the q-Gaussian distributions: 
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arising from the maximization of the Tsallis entropy 



3, 0, Q: 



19j. Nevertheless some 



authors pointed out that this approach is too limited in scope, since q-Gaussian are too 



specific to 



land 



systems [15 
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e the richness of quasi-stationary states observed in long-range interacting 
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22l | . More recently some effort was made to generalize the Central Limit 



Theorem (CLT). It was shown that for a special type of correlation (ca l 



the sum of stochastic variables tend to a q-Gaussian distribution [23 
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ed q-independence) 
25| . Some attempts 



were made to provide specific examples where such type of behavior occur 26 
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33| . In a recent paper Hilhorst and Schehr 3J] have shown analytically 



28|, 



that the sum of correlated stochastic variables for two such examples (see Ref. 271 and 
Reference 12 of [34]) are in fact not q-Gaussians, but functions that are only approximated 
by q-Gaussians. Dauxois also raised some additional points on the applicability of this 
approach 35|,l36||. Rodriguez and collaborators described a simple model with scale invariance 



where this generalized CLT indeed applies [37[. This shows that examples satisfying its 
hipothesis are difficult to obtain. In this work we are interested to study the sum of stochastic 
variables given by the momenta of particles in a classical long-range interacting system, the 
31, 32], and possible relationships with the generalization of the CLT, and 



HMF model 
to ergodicity. 

The HMF model is a paradigmatic and well-studied model, displaying many charac- 
teristics of long-range interacting systems. In references it is argued, based on numerical 
simulations, that the q-generalized CLT is realized in the sum of velocities of particles of 
the HMF model. On this basis it was also argued that the system is non-ergodic. Pluchino 
and co-workers refined this assertion by pointing out that q-Gaussians are only obtained in 
some realizations of their numerical simulations (depending on the microscopic details of the 
initial conditions), while in the remaining cases the distributions are either Gaussian or a 

n 

sort of intermediate function between Gaussians and q-Gaussians '3^^ . In this paper we show 
that if an attractor ever exists in the distribution space it is a Gaussian. The distributions 
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are neither q-Gaussians nor 



obtained in the conditions used in references 
limit distributions in the sense of a genuine CLT. 

We also discuss that the time averages over a single particle history and ensemble averages 
are different, albeit diminishing with time. This leads to the conclusion that even-though 
the system becomes more ergodic (mixed), the quasi-stationary state relaxes to the final 
thermodynamic equilibrium before the ergodicity condition is reached. This explains the 
failure of Lynden-Bell theory of violent relaxation, which requires a complete mixing state. 



II. CENTRAL LIMIT THEOREM AND LONG-RANGE INTERACTIONS - THE 
HMF MODEL 

Let us start by writing down the hamiltonian for the HMF model j^: 

N 2 I N 

^ = E| + ^E[i-co<^^-%)]- (2) 

i=l j,i=i 
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This model describes a system of rotors with angular positions 6i and angular momenta 
(velocities) Pi. The magnetization M is defined as the average value of the modulus of the 
vector (sin 6', cos^^): 

1 ^ 
= ITfUsin^i, 

i=i 

1 ^ 

i=l 



M = ^Ml + Ml- (3) 

The classical central limit theorem states that the sum of n non-correlated stochastic 
variables with finite standard deviation converges to a Gaussian distribution. The standard 
deviation of the summed variable scales with ^/n 38|] . Let us then consider the discrete time 



average momentum of a fixed particle in a system with N particles, say the k-th particle: 

1 " 

where At is a fixed time interval. This average is also a stochastic variable when considering 
its values for each of the particles in the system. If At is sufficiently large so that 
the momenta of the same particle at different times are not correlated, and for statistically 
stationary states, the standard deviation scales as n~^/^. Due to the mean- field nature of the 
system, all relevant physical quantities can be obtained from the one-particle distribution 
function /(r, p, t). If the standard deviation of pj^ approaches zero for large enough n, 
then computing averages using /(r, p, t) is equal to time-averages over one single particle 
evolution, i. e. the system is said to be ergodic. For the purposes of the discussion below and 
for comparison with other works, we also define the reduced average momentum (vanishing 
average and unit standard deviation) by: 

where a is the standard deviation of p^: 



1 " 



.^Ep-I. (6) 



and 

1 



1 ^ 



which vanishes for systems with zero total momentum. A (central) limit exists for the sum 
of stochastic variables if the ensuing distribution has a limit, or rephrasing, if any finite 
statistical moment has a limit for n oo which then corresponds to the moment of the 
limit distribution. This is the meaning of the word limit in the name of such theorems. 
Reduced variables are useful for comparing distribution functions with different averages or 
dispersions, but which are of the same type, e. g. all Gaussian distributions are the same 
with respect to the reduced variable. Let us note that unfortunately contrary to what is 
stated in Ref. [30|, ergodicity is not equivalent to requiring that the distribution of each 



single t/k and the reduced distribution of the momenta are the same. As a counter example, 
consider a system which randomly jumps between two states represented by —1 or +1. The 
state of the system after the i-th jump is labeled by the random variable Zi. The distribution 
of probabilities after a large number n of jumps converge to P(— 1) = P(+l) = 1/2. Since 
the variables Zi are uncorrelated and have finite standard deviations, the sum = Y^i-Li 
approaches a Gaussian distribution as n increases, which is obviously different from the 
distribution for Zi. Clearly this does not imply that the system is non-ergodic, as this is a very 
simple case of an ergodic system (time averages are equal to ensemble averages). Therefore 
arguing that a fixed time reduced distribution of velocities different from the distribution of 
the summed variable i/i as defined in eq. implies non-ergodicity is not correct. We show 
below that for long-range interacting systems the dispersion in p^^, scales much slower than 
n"^/^, even at thermodynamical equilibrium, due to the presence of strong velocity auto- 
correlations. If the life-time of the quasi-stationary state is shorter than the time required 
for a in eq. ([6]) to significantly approach zero, then the assumption of ergodicity used in 
Lynden-Bell theory breaks down. In the present paper we study this issue in numerical 
simulations of the HMF model. 

We first consider the issue of whether a limit distribution exists for the reduced variable 



i/k, and if so if it is a q-Gaussian, as claimed in Refs. [30|, |3l|, l32|, l33[. In order to extend 



the analysis of distributions we perform numeric simulations using the parameter values 



of [3QL 
Ref. 
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32] but improving statistical precision. Let us consider the case studied in 



301 ] with the following parameters: = 100, At = 40 and n = 50, with a transient 
time tj. = 100 to allow the system to reach the quasi-stationary state, with all particles 
at the origin 6i = (magnetization = 1) and the momenta uniformly distributed in 
an interval —po < p < po, with an energy per particle E = 0.69. The simulations were 
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performed using a fourth-order simpletic integrat or |39[| . the same as in [3(|], and repeated 
using a different simpletic integrator for checking [40]. The maximum relative energy error 




admitted was 10^^, and we average over 10,000 realizations (10 times more than in 30|). 
which gives a much better statistical accuracy in the tails of the distribution. At this point 
we note that the use of a good quality random number generator is strongly required to 
avoid unwanted correlations when generating the initial conditions. For a good discussion 



of this point see 
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42l |. The resulting distribution for is shown in the upper pane 



of 



figure [H together with the best fitted q-Gaussian with q = 1.67 (the value obtained in [30|| is 
q = 1.65). It is clear from the tails of the distribution that it is not a q-Gaussian. The lower 
panel in figure [1] shows the reduced velocity distribution function at time t = 2100 (the final 
time of the interval over which the sum (ISD runs), which coincides with the one obtained 
in 



301 ]. This distribution is not Gaussian and can be well approximated as a limit case of a 
Lynden-Bell distribution for a two- level initial distribution for A4 ^ 1 (see eq. below). 
It is straightforward to show that the distribution of t/k cannot be a q-Gaussian from the 
fact that for g > 1 the statistical moments of the reduced variables yi: 

1 ^ 

M, = ^T.\y>^\', (8) 

k=i 

diverge for j > j' for some j'. In fact eq. ([H]) is an approximation from the simulations for 
the true moments obtained for N oo. If the true value of Mj diverges, then computing 
the simulation value of Mj from eq. ([H]) for increasing would show no convergence to a 
finite value. Nevertheless in all our simulations the value of Mj always converged to a finite 
value for a given j such that the corresponding fitted q-Gaussian implied an infinite value. 

Another important point is that the situation at the final time is not a limit state, i. e. the 
moments of the reduced variable yk do not converge to a constant value for increasing time, 
as seen in fig. [2l Along these lines, for the situation considered in [3^ no limit distribution 
is obtained. 

Extending our analysis let us consider the case studied in ref. jsil, with A^ = 100,000, 
At = 100, n = 10, 000 and E = 0.69, with 50 realizations and the same type of initial 
conditions. The results are shown if fig. [3l together with the q-Gaussian fit with smallest 

error. It is evident from the figure that no q-Gaussian is realized. The lower panel in 
fig. [H] shows the moments M3 and M4 for 15 different realizations as function of time. The 
reference values of the moments for a reduce variable having a Gaussian distribution are 



FIG. 1: Left panel: histograms for the distribution of the variable yi over 10,000 realizations, 
N = 100, At = iO, E = 0.69 and the best fitted q-Gaussian with q = 1.67 and (3 = 5.88. Right 
panel: momentum distril 



FIG. 2: Moments M3 and M4 of the distribution of the variable yk, for the case in fig. [TJ as a 
function of t = nAt (varying n). Here it is clear that no limit reduced distribution is attained as 
must be the case for a limit theorem (see comment after eq. ([7])). 

M3 = 1.5958 and M4 = 3. The figure shows that all cases, except one which is far from a 
limit situation, are converging to a Gaussian distribution. The present case is an example 
of sluggish convergence of a stochastic process to a Gaussian, as discussed in a different 
context in References 



43 



4J|. The convergence rate itself depends on the details of the 



initial conditions. 

From the results presented above we can draw the following conclusions: the only true 
attractor observed in distribution space for the sum in eq. ([5]) is the Gaussian distribution. 



In all cases considered here and in 



251] the distribution is never a q-Gaussian. The 



only situation where a limit exists for the distribution of the sum variable is such that it 
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FIG. 3: Left panel: distributions functions for yi in eq. ^ with N = 100, 000, At = 100 n = 10, 000 
and E = 0.69. Some distributions are very close to a Gaussian. The best fit using a q-Gaussian 
(distribution with minimal error) is also shown. Right panel: moments M3 (lower set of lines) 
and M4 (upper set) as a function of time t = nAt for 15 different realizations showing a sluggish 
convergence to the values for a Gaussian. 

converges to a Gaussian, as show by the time evolution of moments. The claim that a 
generalized central limit theorem leading to a Tsallis distribution is observed for the HMF 
model is therefore invalidated. This do not precludes the applicability of such theorem in 
other situations. 



III. NON-ERGODIC BEHAVIOR 



Now we turn to the ergodicity problem for this system. As explained above, we look 
at the time evolution of the dispersion a in eq. ([6]). We consider the HMF model with 

= 10, 000 particles. At = 40, with a uniform initial distribution in an interval given by 
f{p,9,t = 0) = /o = 1/2^06*0, if —po < P < Po; < 6^ < 6^0, and zero otherwise. The 
parameters po and ^0 can be used to adjust the energy and initial magnetization. The 



45|. In that 



predicted final state being computed using the procedure described in Ref. 
direction we first note that the value A^^. in eq. in the quasi- stationary state can be set 
to zero by a global shift in the angles. The distribution function of the quasi-stationary state 
predicted by Lynden-Bell theory for an initial distribution with two levels (/o and zero) is 



given by 



13|: 



fLB{p,e) 



1 _|_ Q-|3e{pfi)-^l ■' 



(9) 
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where 

2 

e{p,9) = ^-MySm9, (10) 

is the one particle energy, and (3 and /i are Lagrange multiphers determined from the energy 
and normahzation constraints: 

J e{p, e)dpde = u, (11) 

/ fLB{p,0)dpde = l, (12) 
together with the corresponding expression for Al^: 

My = J hB{p, 0) sine dpde. (13) 

Equations fill I) and 0121] form an algebraic system in the unknowns (3 and fi and can be solved 
numerically to explicitly determine Jlb- Is is thoroughly discussed in the literature that 
the predictions of Lynden-Bell theory yield at best approximate results due to incomplete 
mixing. We will show that this is in fact due to the strong velocity auto-correlation which 
prevents the system to reach an ergodic state, in the sense discussed above. We show in 
fig. m the dispersion (standard deviation) for the summed variable p^, for some values of the 
magnetization of the initial state for = 10, 000. The life-time of the quasi- stationary state 
can be inferred from fig. [5] showing the kinetic and potential energies per particle for the case 
A4 = 0.6 (the other cases has essentially the same life-time) as r ^ 2 x 10^. The t^^/^ scaling 
for a is reached only after a very long time, of the order of the life-time of the stationary 
state. The system evolves to the thermodynamical equilibrium well before the dispersion 
a can attain a significantly small value. Figure El shows the velocity distribution function 
for two times compared with the Lynden-Bell distribution function for Ai = 0.6. Even 
though the distribution is slowly approaching the theoretical prediction while a decreases, 
the difference remains always significant. 

IV. CONCLUDING REMARKS 

We have shown that the HMF model is non-ergodic due to strong auto-correlations in 
particle velocities, which are also responsible for anomalous diffusion and Levi walks, re- 
straining the dispersion in the time average of the velocity of each particle to approach 
zero before the system finally settles down in its thermodynamical equilibrium. The system 
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FIG. 4: Dispersion a of the variable Pf. in eq. ^ for different initial magnetizations. The system 
first evolves 1,000 units of time, allowing the the violent relaxation to occur and settle down close 
to a stationary state. The time indicated in the figure is measured from the end of this transient. 
The dashed line ^ 
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FIG. 5: Kinetic (K) and potential (V) energies per particle for the case = 0.6. The transition 
from the quasi-stationary state into the final equilibrium (Gaussian) state is given by r = 2 x 10^. 

has not enough time to reach ergodicity and therefore any intermediate state before final 
equihbrium is incompletely relaxed. Ergodicity is only fulfilled for asymptotic times, well 
after thermodynamic equilibrium occurs. This clearly points out that a solution of the vi- 
olent relation problem must necessarily rely on a dynamical theory, based for instance in 
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FIG. 6: Distribution functions from simulation for the case A4 = 0.6 at two different times com- 
pared to the prediction of Lynden-Bell theory. The Lagrange multiphers obtained by solving 
eqs. (dl]) and are /3 = 4.86 and /i = -6.67. 

the solution of a proper kinetic equation as already pointed out by Chavanis [l^. Many of 
the points presented here can be extended to other systems with long-range interactions as 
long as strong velocity auto-correlations are present. We have also shown, if simulations are 
performed with care to statistical precision, that the only attractor observed for the sum of 
the velocity of a single particle at fixed time intervals is the Gaussian distribution, thus in- 
validating the claims by Tsallis and collaborators that such attractors should be q-Gaussian 
distributions. 
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